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1. INTRODUCTION 



Unsteady flow in the passages of wave rotor devices can adequately be 
modelled on a one-dimensional basis. However, this modelling can be quite 
involved due to the peculiar characteristics typical of wave rotor type flows. 
The numerical calculation has to provide approximate solutions of 
time-dependent compressible fluid flow problems which involve discontinuities 
and strong wave interactions. Ref. (1) lists three criteria which such 
approximate solutions should satisfy simultaneously: (i) the solution must be 

reasonably accurate in smooth regions of the flow. Continuous waves 
(rarefaction waves, compression waves) should propagate at the correct speed 
and should maintain the correct shape which involves steepening or spreading 
at the correct rate; (ii) discontinuities which are transported along 
characteristics (gradient discontinuities, contact surfaces), should be 
modelled by sharp and discrete jumps, and should be transported at the correct 
speed; and (iii) nonlinear discontinuities such as shocks should be computed 
stably and accurately. 

In addition, the complex pattern of shock waves and contact surfaces that 
could evolve in wave rotor devices precludes the use of numerical metluods 
wiiich rely on either some type of artificial viscosity or a special tr:MCment 
of discontinuities. Such methods would quickly h>ecome quite impractical for 
this application due to programming difficulties and cost of execution. 

Computation of such solutions has generally bt^en carried out by solving a 
set of finite difference equations which approximate the governing 
differential equations of flow. All such schemes inherently have a finite 
amount of dissipation as well as dispersion of the wave modes they model, and 
it is difficult to construct difference schemes which simultaneously satisfy 
the criteria given above. Stability problems may also be an added concern for 
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these schemes 



In view of the foregoing, an alternative approach to solving wave rotor 
type flows was sought, and the purpose of this report is to describe such a 
scheme along with some results. The scheme is known variously as Gliram*s 
method, the Random Choice Method (RCM) or the piecewise sampling method. The 
method evolved from a constructive proof of the existence of solutions to 
systems of nonlinear hyperbolic conservation laws given by Glimm (Ref. 2). 
Chorin (Refs. 3 and 4) developed the scheme into an effective numerical tool 
for gas dynamic applications, with emphasis on detonation combustion problems 
and reacting gas flows. Although the ROM computes solutions on a fixed grid, 
it is not a difference scheme, utilizing solutions of locally defined Riemann 
problems as the basic building blocks for the global solution. Each of the 
local Riemann problems (defined in more detail in section 2) provides an 
analytically exact elementary similarity solution. By means of a suitable 
sampling procedure, usually of a pseudo-random or quasi-raniiom nature, the 
similarity solutions are superposed to construct the approximate solution to 
the equations. 

With an appropriate sampling technique, the RCM in one dimension is 
possibly superior to any finite difference scheme in meeting the criteria 
established above. 
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2* METHOD 



2, 1 . Solution Procedure 

The method models the one-dimensional, compressible, inviscid Euler 

equations, expressed in conservation form as 

au 3F(U) , 

at 3x 



U(x,t) = 


P 

pu 


and 


F(U) = 


^P 

9 

PU“ -H p 




^E ^ 






(E + p),. 



Here E is the total energy per unit volume and may be expressed as (for a 
polytropic gas) 



E 







e ^ internal energy per unit mass 

= _L/p\ 



Y-l j 

p is the density, p is pressure and u is velocity in the one space 
dimension being considered here. With initial data specified in the form 

U(x,0) = 9 (x) , 

an initial value problem is defined for the Euler equations. The simplest 
initial value problem for which discontinuities appear is the Riemann problem 
to find the gas flow resulting from an initial state in which the gas on the 
right of an ’origin* is in a constant state, and the gas on the left is in 
another constant state, i.e., 



with 



9(x) 



Ul 

Ur 



X < 0 
X > 0 



Ul.r = 



►'L,R 

(|Ju)l,R 

*^L,R 
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where the subsripts L and R denote the left and right sides of the 
‘origin*, here arbitrarily prescibed at 0 . That is, the Riemann problem 

consists of prescribing constant initial data on either side of an origin 
where a jump discontinuity exists. As mentioned before, the solution of the 
problem constitutes a basic building block of the random choice method. A 
special case of the Rieraann problem in which ul = u[^ = 0 is often referred 
to as the shock tube problem. The answer to the problem is that there are 
four possible types of subsequent flow, depending on the inequalities in the 
left and right side data prescribed. Thus, in both directions from the 
origin, a shock or a centered rarefaction wave may propagate, giving rise to 
the above mentioned four different possibilities. Fig. (1) illustrates the 
special case of shock tube type flow and the evolution of the wave pattern. 

Fig. (2) shows the simple fixed Cartesian grid set up for the method. 

Let Lx be a spatial increment and At a time increment. The solution is to 

be evaluated at time (n + l)At , n being a non-negative Integer, at 
spatial increments iAx , i = 1,2,3, . . . The initial data is prescribed 
for each time step at nAt in a piecewise constant manner i.e., it consists 
of intervals of length Ax where the data is constant, separated by jump 
discontinuities : 

U(x , nAt) = , (i-“)Ax < x < (i-hy)Ax 

The solution at time (n+l)At then is required to have the same property, 

i.e., it is piecewise constant over an interval Ax , and it serves as the 
initial data for the next time step: 

U(x,(n+l)At) = , (i - -j)Ax < x < (i-K^)Ax 
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This procedure defines a sequence of local Riemann problems to be solved at 
each time level. On the grid shown in Fig. 2, for example, initial data would 
be specified ai: points I, 3, 5 setting up a succession of Rieraaun 

problems defined by each pair of states (1,3) , (3,5) , (5,7), with the 
discontinuities at the midpoint of each, i.e. , at 2, 4, 6, .... etc. If the 
time step increment At is calculated such that 

At < o.(Ax). max (|u?|+a^) , with 

0 < o < i 

then the waves generated at the discontinuities of adjacent Riemann problems 
will not interact, as shown schematically in Fig. 2. 

Each of the local Riemann problems yields an exact analytical solution, 
with the resulting wave structure a particular corahinat ion/ v^ar iat ion of the 
general structure shown in Fig. 3. 

In the x-t plane, the solution to a Riemann problem consists of 
essentially four regions connected by three waves. Thus states I and TV are 
the prescribed left and right states for the problem, and states II and III 
are the ’starred* middle states separated by a slip line or contact 

d X 

discontinuity = u* . The velocity, u , and pressure, p , are 

continuous across the contact, but in general is not. Thus Uj^.* " > 

PL* = PR* PL* ^ ^R* . S[,b , S2,b ^L>f * ^2>f represent 

respectively the backward and forward facing waves generated at the point of 

discontinuity and may be either shocks or rarefaction waves. 

Still referring to Fig. 3, it is seen that at a time nAt < t < (n-M)At, 
the exact solution of the local Riemann problem for the interval ((i-l)Ax , 
iAx] may actually consist of several distinct states. Consider now a 
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translation of each interval [(i-l)ax, iAx] to ^ J such that 

the discontinuity (i.e., the point from which the waves are generated) is 

centered at a zero origin. Let 0 be the value of a random variable » defined 
over the interval ^ 



C = OAx , i.e. - < 5 < + ^ 

Also, define (x,t) , aAt < t < (n+l)At , to be the exact solution to 

exact 

each Riemann problem. Using the value of ? to fix a point in the interval 
Ax of each Riemann problem, the exact solution at that point is determined 
and assigned to either the left or the right grid point, depending on whether 
^ is < or > 0 . Thus, if the point fixed by C is F* (Fig. 3), the 
exact solution to the Riemann problem at that sampled location is assigned to 
the grid point on the right and if the sampled point is P'* , the solution at 

that location is assigned to ihe grid point on the left, i.e., for a typical 
interval [(i-l)Ax, iAx] , 



if C < 0 , U"'*’} = ^ (^ , t) 

i~l exact 

and if C > 0 , ^ (^ , t) 

1 exact 

It is seen immediately that although the solutions are computed on a 
grid in this method, it is not a differencing scheme. Also, instead of using 
a weighted average of the Riemann problem solution to arrive at the solution 
for a grid point t, the RCM samples a particular value from an explicit wave 



t The Godunov method, for example implements 




^ i^(i+j)Ax 

Ax / , , I . . 

(l-j)AX 



,n-F 



exact 



(x , t )dx 
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solution, thus eliminating the smoothing out of wave transport and interaction 
information inherent in averaging. This leads to the ‘infinite* resolution of 
contact discontinuities and shocks that the scheme displays. 

From the foregoing discussion, it is evident that the success of the 
scheme hinges, to a large extent, on the inexpensive and exact solution of 
Rieraann problems and an appropriate sampling technique. Ref. (3) describes a 
modification to an iterative method due to Godunov (Ref. 5). TheoreticaL 
details for the Riemann problem solution are also given in Ref. (6). 

The mathematical properties required in a sampling procedure applicable 
to this scheme are defined in Ref. (1). A brief description of the procedure 
is given below. 

In previous computations using the RCM, random sampling with some 
variance reduction technique (stratified sampling); was ust-d, i.e., tlie values 
wer<i taken from the random number generator installed in the computer (Ref. 

3). It was shown in Ref. (1) that a more accurate form of sampling is a 
technique due to van der Corput (Ref. 7). The sequence gener<ited is, strictly 
speaking, non-random, but has particular statistical properties that are 
suitable to the application. The sequence is referred to as quasirandom and 
is generated as follows: 

The binary expansion of natural numbers may be expressed as (witli R=2): 

n = AqRO + A^rI + A2r2 + + Aj^Rm , (0 ^ A [. < R) 

m 

i.e. n = I A^^.2^ , with = 0 or 1 , n = 1, 2, 3, ... 

k=0 
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Next, the digits of the binary numbers are reversed and a decimal point is put 
preceding the number; this gives the numbers 

4-n = AqR'^ + AiR"2 + ... + 
m 

or, ^ Aiq. , again with A^ = 0 or 1 

k=0 



Conversion to the decimal scale of these numbers yields the required sequence 

of quasirandom numbers defined over the interval [0,1], i.e., 

if-n (decimal) = Gn + y 

or Of! = (decimal) - 



and ” Oj^.Ax as defined earlier. 

The first few elements of the sequence given below illustrate the construction 



11=1 (decimal) = 
2 

3 

4 

5 

6 

7 

8 



I (binary); 
10 

II 
100 
101 
no 
111 
1000 



= 0.1 (binary ) 

0.01 
0.11 
0.001 
0. 101 
0.01 1 
0.111 
0.0001 



= 0.5 (decimal) 

0.25 
0.75 
0.125 
0.625 
0.375 
0.875 
0.0625 



The van der Corput sequence is 'equidistributed ' , and yields better results 
than those obtained using a 'stratified' random sampling technique. 

The subroutine employed in the program to compute the random numbers is 
described in Appendix B. 

2.2. Boundary Conditions 

In general, the implementation of boundary conditions in the RCM is 
quite straightforward, but does require some thought. Referring to Fig. 2, 
the h.c.'s are specified at points 1 and N for the left and right boundary 
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respectively. Note that if the sampled solution at (n-K)at corresponds to a 
random number ^ < 0 , the solution is assigned to the grid point on the 
left. For the Riemann problem defined by points 1 and 3, the sampled solution 
would then be assigned to grid point 1 at (n+l)At ; however this is 
overridden by assigning the proper boundary condition at 1 again, and there is 
no contradiction. A similar procedure is adopted at the right hand boundary 
when Cn > 0 • 

The subroutines for the boundary conditions are named in the format 
BCxn , BC standing for ^Boundary Condition, x being either L (for _Left), 
or R (for _Right boundary) and n being a number from I to 5 with the 
following designations: 

1 - solid wall condition 

2 “ outflow at constant static pressure 

3 - special formulation (* piston* inflow) 

4 - isentropic inflow from reservoir 

5 “ special formulation (rarefaction wave cancellation) 

2.2.1. Solid Wall Conditions 

The solid wall boundary condition requires a zero normal velocity at 

the wall for inviscid flow computations. Due to the random sampling involved 
in the method and the lateral movement of the sampled solution to the 

left or right of the discontinuity, the condition is difficult to implement 

uniquely. However, the procedure adopted here is found to yield reasonably 

accurate results for the applications intended. (Note that the difficulty is 

not unique to this method only. The implementation of zero mass flux through 

a surface is difficult per se for the Euler equations). 

Referring to Fig. 2, let the physical boundaries be at point 2 and 
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point (N-1) for the left and right sides respectively. However, the boundary 
conditions are specified at point I (point N) for the left (right) side as a 
fictitious *mirror * state of the conditions at point 3 (point (n-2)) 
respectively, but with the reverse sign taken for the velocity component. 

Thus, for the left hand boundary Rieraann problem, 

= P(3) , PL = P(3) , UL = -u(3) 

PR = P(3) , PR = p(3) , UR = u(3) 
and, analogously, for the right hand boundary Riemann problem, 

PL = p(N-2) , PL = p(N-2) , UL = u(N-2) 

PR = p(N-2) , PR = p(N-2) , UR = -u(N-2) 

The solutions are then sampled in the manner outlined earlier. 

2.2.2, Outflow Conditions 

For subsonic outflow, only the static pressure p is defineci, with 
the continuation condition being applied to the rest of the variables. Thus, 
for the right hand boundary for example, the Riemann problem is defined as 
follows : 

PL = p(N-2) , PL = p(N-2) , UL = u(N-2) 

PR = Pout > PR = p(N-2) , UR = -u(N-2) 
where Pout the specified outlet pressure. If the flow going out is 

supersonic, there can be no propagation of disturbances upstream, and the 
continuation condition is implemented for all the variables, i.e,, the Riemann 
problem now is the trivial case defined by 

PL = p(N-2) , PL = p(N-2) , UL = u(N-2) 

PR = p(N-2) , PR = p(N-2) , UR = -u(N-2) 

2.2.3. Special Formulation of *Piston* Inflow 

In general, for idealized wave rotor flows, hot combustion gases are 
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introduced into the rotor through nozzles angled such as to allow the flow to 
*slip onto* the rotor, i.e. , without incurring incidence or deviation angle 
losses. Also, in the ideal treatment, the air in the passages of a wave rotor 
is exposed to the hot gas at high pressure instantaneously. The idealizations 
allow for uniform conditions to be prescribed at the hot gas inlet port. 

Thus, a *special* form of inflow boundary condition needs to be specified 
here, namely, the static pressure, the velocity and the density of the 
incoming hot gas. Although equivalent to specifying the total pressure and 
temperature in the usual inflow boundary condition treatment, some thought is 
required in wave rotor type flows when specifying Pgas > ^^gas '^gas 

This is because only a shock wave needs to be generated, with no waves 
travelling opposite to the direction of flow. The solution to the Riemana 
problem would then consist of just two states connected by a single shock 
wave. The flow is equivalent to that generated when a piston is pushed 
instantaneously into a gas at rest. In general, the state of the air inside 
the rotor passage is known; explicit relations for two states connected 

through a shock wave are given in Ref. (6). These so-called transition 

functions help in specifying the boundary conditions for the incoming flow 
properly. 

If we consider the left boundary for this inflow, the Riemann problem 
is set up as: 

PL “ Phot gas ’ ^-'L “ Phot gas > " ^hot gas 

Pr = p(3) , PR = p(3) , UR = u(3) 

with pl , P]_, and u^ having been chosen in accordance with the 



considerations discussed above. 



2. 2,4. Isentroplc Inflow From Reservoir 

The induction of fresh charge or air onto the rotor usually 
corresponds to an isentropic inflow situation. The flow in the vicinity of 
the passage end can be treated as quasi~steady , with the assumption that no 
flow separation takes place when the flow enters. Two boundary conditions are 
required for this type of inflow; these are provided by the conservation of 
energy in the flow from the external region to the inlet (assumed to be 
steady), and by the prescibed entropy level of the gas in the external 
region. 



where the subscripts ’in’ and ’tot’ apply to conditions at the inlet of the 
passage and external reservoir respectively. The sonic velocity is denoted by 
a , and flow velocity by u . Note that knowledge of the Riemann variable 
arriving at the passage end from within the passage is required to be able to 
solve the energy equation above for a^j^ and u^i^ . For the left end, for 
example , 



The boundary conditions may thus be expressed as 




2 2 

r a 

Y-1 tot 



^in ” ^tot 




which together with the energy equations yields 




Y+1 



Y-1 



and subsequently the other variables 
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The simple analytical treatment given above has to be modified 
somewhat if a contact discontinuity is formed when the inflow begins. This 
is due to the fact that the value of the arriving Rieraann variable is changed 
across such a discontinuity, which thus leads to an additional unknown. 
Procedures for solving the inflow for these situations are given in Ref. (8). 
In the program developed here, reasonably good results are obtained by setting 
the velocity at the boundary point equal to the velocity at the point nearest 
the physical boundary. For the left end e.g., the variables for the left 
state of the Riemann problem are obtained as follows: 

u(l) = u(3) , 

a reasonably accurate assumption just at the point of inlet opening. 

Then, from the ’energy ellipse’. 



a(l) 





u(l)2 



M(l) 




incoming Mach number 



p(l) = 



-Etot, 



[I 



Y-i 



M(l)^ 



y/(y-i) 



with similar isentropic relations to compute other flow variables. Note 
that once the interface or contact discontinuity has moved i\ certain distance 
inside the passage, the simple analytical expressions given earlier in the 
section can be used, since now the value of the arriving Riemann variable 
would be known at the boundary. 

2.2.5. Special Formulation for Rarefaction Wave Cancellation 

The spreading of rarefaction fans leads to unwanted wave reflections 
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which occupy large zones in the passages of wave rotors • Fig. (4) shows a 
wave diagram proposed by Spectra Technology, Inc., which incorporates 
so-called *wave management* or ’tuning* ports to ideally cancel (and otherwise 
attenuate) impinging rarefaction fans. The physical boundary conditions are 
thus dictated by the flow developing in the passage, i.e., the port has 
non-uniform flow conditions in it, which at each point match those of the flow 
at the end of the passage so as to disallow any reflections to take place. 
Numerically, this condition is achieved by implementing the continuity 
condition across the boundary for all the flow variables involved. For the 
left boundary, thus, the Riemann problem is defined by: 

PL = P(3) , PL = P(3) , UL = u(3) 

PR = P(3) , PR = p(3) , UR = u(3) 

and analogously for the right boundary. Note that these boundary conditions 
may involve either inflow or outflow. 

2. 3 Example Calculations 

The listing of the program is included in Appendix A, and the various 
names for the variables are listed in Appendix B, along with some instructions 
on how to use the program. No effort as yet has been made to optimize the 
code either for storage requirements or for execution efficiency. 

In this section, some sample calculations are carried out using the code, 
to illustrate its usefulness in constructing idealized design point wave 
diagrams which can serve as the starting configuration for detailed 
construction of diagrams incorporating real flow effects. 

2.3.1. Test Case for 1-D, Inviscid, Unsteady, Compressible Flow 

Fig. (1) illustrates the initial conditions in a shock tube, with the 
diaphragm at xq • Sod (Ref. 9) suggested a test case for hyperbolic 
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conservation laws with the following conditions as initial states in the shock 



tube : 

PI = KO , PI = 1.0 , u[ = 0.0 

P5 = 0.1 , P5 = 0.125 , u5 = 0.0 

i.e., the gas on either side of the diaphragm is in a quiescent state 
initially. The ratio of specific heats is chosen to be 7/5, and Ax is 
chosen to be 0.01. 

The solution (before any wave has reached either the left or right 
end) is shown in Fig. (5). The squares shown at locations , x2 , X3 

and X 4 in the density plot give the analytically calculatc^d amplitude and 
location of the head - and tail waves of the left-running r.iref act ion , the 
contact surface moving to the right and the shock wave moving at supersonic 
velocity to the right respectively. The solid lines are the solutions 
obtained by the RCM at different time levels; the zero numerical diffusion 
feature of the method is evident in the ’infinite* resolution of the contact 
discontinuity and the shock, and the dispersion (phase error) is within one 
grid spacing. The constant states are perfectly realized. 

It is these features of the method that make it very attractive for 
application to wave rotor type flows, since the successful design of the 
device is predicated on being able to accurately compute wave arrival times at 
the various ports. 

2.3.2. Wave Turbine Experiment 

Ref. (10) describes the wave rotor experimental set up at the 
Turbopropulsion Laboratory. Initial tests being carried out currently are 
with the wave rotor in a turbine mode, i.e., one side of the rotor is blocked 
off, and high pressure air is brought onto the rotor and taken off again from 
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the other side. The passages of the rotor being angled at 60® to the axis, 
the 180® reversal in the direction of the fluid flow creates an angular 
momentum change, in turn generating large turbomachinery work coefficients. 
Fig. (6) shows the wave diagram computed using the code. The movement of the 
rotor is from top to bottom. At t=0 , the high pressure air is brought 

into contact with quiescent atmospheric air in the rotor passages, at point a. 
This corresponds to the ’piston’ inflow boundary condition described in 
section 2.2.3.. A shock, S , is generated immediately, (idealized case of 
instantaneous cell opening), which travels from the right to the left, and 
strikes the solid wall at the left end. The reflection of the shock takes 
place at point b according to the solid wall boundary condition described in 
section 2.2.1.. Behind this shock, and moving at a slower velocity is the 
contact surface, I , which penetrates into the passage only a fractioi^al 
distance before encountering the reflected shock, RS , at point c. The 
reflected shock is transmitted through the contact surface, (bringing the flow 
to a near halt), and reaches the right side at point d, whereupon the inlet 
port is closed. The air trapped in the rotor passages is now at a high 
pressure and in a quiescent state. When this air is released at point e to a 
low pressure region, a rarefaction wave is generated, R , which travels to the 
left, spreading out in the process. It interacts with the stationary contact 
surface, I , setting it into motion again, and reflects off the solid wall at 
the left as RR . The boundary condition imposed at point e is the outflow 
at constant static pressure condition described in section 2.2.2.. The outlet 
port is closed at a time when the exit velocity falls to about half its 
initial value. 

This experiment embodies two fundamental processes in wave rotors: 
those of cell filling and cell emptying. Almost all the other processes 
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typical to wave rotors are combinations of the cell filling and cell emptying 
unit processes. Comparison of the ideal computed numbers obtained here with 
experimental data will provide information helpful in the identification and 
sources of losses. 

The program is set up to start at t=0 in this case, with initial 
data provided along the entire passage, i.e., from x=0 to x=0. 1863m (the 
actual length of the wave rotor being tested). Since the passages have 
quiescent atmospheric air in them at t=0 , the initial data, of course, 

describes these conditions. Switches for the left and right boundaries 
describe what type of boundary conditions prevail and direct the program to 
the appropriate subroutines. These switches, designated SWL and SWR , for 
left and right respectively, are assigned integer number values which 
correspond to the numeric value of the particular boundary condition they 
represent. Thus, if the left boundary is a solid wall, SWL=1 , 
corresponding to the boundary condition subroutine BCLl . In this example 
then, the initial switch settings at t=0 are SWL=1 and SWR=3 , 
corresponding to a solid wall at the left and a ’piston* inflow at Che right 
(which starts at t=0 at point a). At point d , the switches are reset to 
SWL-l and SWR=1 due Co the closure of Che inlet port. At point e , the 
switches are SWL=1 , SWR=2 , signifying opening of the exhaust port with 
outflow at a constant static pressure. The whole wave diagram can thus be 
packaged into a ’module* subroutine and called from the main program with a 
single call statement. This type of modularity allows for wave diagrams of 
different ’families’ to be developed by simply calling the right ’module’ 
subroutine. 

The next two examples illustrate this concept as they deal with two 
very different types of wave diagrams. 
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2,3.3. General Electric Wave Engine 



Fig, (7) shows a schematic of the wave diagram constructed for the 
G.E. wave engine. Briefly, the device is configured for a gas generator mode 
of operation, with counterflow scavenging, and is capable of producing net 
shaft power. For a fuller description of the machine, see Ref. (11). In this 
example, fresh charge (air) is induced into the rotor (from an external 
reservoir) through the wave action of the rarefaction fan originating at the 
exhaust port opening. The usefulness of the rotor is gauged by the net 
pressure rise across the machine, i.e., the ratio of the total exhaust 
pressure to total (fresh air) inlet pressure. 

For performance estimation purposes, it is sufficient to investigate 
only the exhaust and induction processes as shown in Fig. (8). The initial 
data specified is as follows: the exhausting pressure ratio Pe/Po > 

total pressure ratio across the rotor Pte/Pta assumed total 

temperature ratio this particular cycle, the amount of fresh 

charge induced in is ideally equal to the gases exhausted out, i.e., 

^in ~ ^ut » this mass balance is carried out after each computation to 

correct the assumed temperature ratio I'te/'^ta (which otherwise constitutes 
overspecification of the initial conditions). 

The calculation starts at t=0 , with initial data consistent with 

the chosen pressure and temperature ratios specified along the passage length. 
Initial switch settings are SWL=1 and SV\/R=2 for the solid wall boundary at 
the left and the exhaust to a constant pressure at the right. As shown in the 
figure, a rarefaction fan is generated, propagating to the left and reflecting 
off the solid wall. At time t=ii ^ ^he pressure at the wall has been 
reduced to that outside the passage, pta > which is when the inlet port is 
opened. The switches are now set to SWL=4 and SWR=2 for isentropic inflow 
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from an external reservoir at the left, and still outflow at a constant 
pressure at the right. The exhaust port is closed at time t=i2 which 
corresponds to the exit velocity having dropped off to approximately half its 
steady state value at the beginning of the exhaust process. Now the switches 
are set to SWL=4 and SWR=l , for the solid wall condition at the right. 

The sudden closure of the exhaust port generates a ’hammer^ shock travelling 
to the left, interacting with the incoming interface (shown by dashed line), 
and reaching the passage end at t=i 4 at which time the inlet port is closed, 
with the switches being reset to SWL=1 and SWR=1 . Note the reflected 
shock travelling from left to right generated at the interaction of the 
contact surface and the hammer shock. 

Once this solution is obtained, integration of the mass flux through 
the inlet and exhaust ports is carried out and if the two numbers do not 
match, the assumed temperature ratio is adjusted in the initial 

data, till such time as m^^ = nioyt 

This calculation is sufficient for performance analyses: if the 

entire wave diagram has to be worked out, then at a time t > Cj , hot gas 
from the combustion chamber is brought onto the rotor (the boundary condition 
corresponding to ’piston’ inflow) on the right hand side. This would 
generate the shock to compress the induced air and when this shock reached 
the left end, the transfer port (see Fig. 7) would be opened for such time it 
takes for the compressed air to be completely scavenged out of the rotor. 

Fig. (9) shows some performance curves obtained using the procedure outlined 
above. In Figs. (10a, b, c) are shown three sets of flow parameters at 
different time steps corresponding to the inlet port just opening, the 
exhaust port closing and the inlet port closing; the qualitative distributions 
of the flow parameters in the passage are immediately seen to be accurate when 
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compared with the wave diagram shown in Fig. (8). Of interest is the set of 
plots for the time step when the inlet port has just been closed. The flow 
between the end of the passage and the location of the interface is seen to 
be quite non-uniform in the density plot. At the same time, the shock 
reflected from the interface has reached the right side and reflected off the 
solid wall. These considerations help to decide optimum port opening and 
closing times. For example. Fig. (11) shows what happens if the inlet port is 
not closed at just the time the shock reaches the end, but rather at some 
short time later. The shock now sees an open boundary and reflects off as an 
expansion to match the high pressure behind it with the incoming total 
pressure which is at a lower value. This reflected expansion is manifested in 
the pressure, density and velocity plots of the figure. 

The Entire sequence of wave interactions of this example is computed 
by the RCM without the implementation of artificial viscosity or artificial 
compression methods, or tracking and capturing schemes. This ’hands off’ 
feature of the method renders it eminently useful for fast preliminary 
evaluations of complex wave diagrams for the application at hand. 

The next example computes an idealized wave diagram for the nine-port 
pressure exchanger concept proposed by Spectra Technology, Ref. (12). 

2.3.4. Spectra Technology Pressure Exchanger 

Fig. (4) shows the ideal wave diagram for the nine-port pressure 
exchanger. This configuration is a good case example to compute with the RCM 
because of the different types of boundary conditions that need to be dealt 
with in the evaluation of the cycle. The computation is started at t=0 , at 

the point of high pressure hot gas inlet (driver gas inlet). In the manner 
described in the G.E. wave engine example, the initial data is prescribed for 
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the entire passage at this time step and the boundary condition switches are 
initially set at SWL=1 and SWR=3 for the solid wall at the left, and the 
'piston* inflow at the right hand end. Since there is a multiplicity of types 
of boundary conditions, e.g., three outflow ports, an index, JCOUNT, is used 
to ensure proper sequencing of the switches. The following table is 
presented as an example of the settings of the switches to carry out 
calculations for one cycle. The inflow and outflow port conditions are those 



proposed by Spectra Technology for their idealized diagram. 



TIME STEP, N 


JCOUNT 


SWL 


SWR 


REMARKS 


0 


0 


1 


3 


CYCLE STARTS. HP GAS INLET PORT OPENS 


500 


1 


2 


3 


HP AIR OUTLET PORT OPENS 


1408 


2 


2 


1 


HP GAS IM.ET PORT CLOSES 


1765 - 


3 


5 


1 


HP AIR OUTLET PORT CLOSES. TUNING PORT 
LI OPENS 


1816 


4 


2 


1 


TUNING PORT LI CLOSES. IP GAS OUTLET 
PORT OPENS (PORT El) 


2069 


5 


2 


5 


TUNING PORT R1 OPENS 


2261 


6 


2 


1 


TUNING PORT R1 CLOSES 


2595 


7 


5 


1 


IP GAS OUTLET PORT CLuSES. TUNING PORT 
L2 OPENS 


2636 


8 


2 


1 


TUNING PORT L2 CLOSES. LP GAS OUTLET 
PORT OPENS (PORT E2) 


3029 


9 


2 


5 


TUNING PORT R2 OPENS 


3237 


10 


2 


4 


TUNING PORT R2 CLOSES. LP AIR INLET 
PORT OPENS 


4961 


11 


1 


4 


LP GAS OUTLET PORT CLOSES 


5529,0 


0 


1 


3 


LP AIR INLET PORT CLOSES. CYCLE 
COMPLETED 



The total cycle time as calculated by the RCM is 3.0676 mseconds, which 



compares well with the time computed by Spectra Technology (using the 
FCT-SHASTA algorithm) of 3.07 mseconds. The execution time on an IBM 
370-3033AP for the 5529 steps computed in the example above was 3 minutes 
38 seconds, including the I/O operations and the graphics. 
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Figs. (12a, b, c) show three sets of plots of the flow parameters for 
the following cases: a) the H.P. air outlet port opens on time, i.e., just as 

the shock reaches the left end of the passage, b) the port opens before the 

shock has reached the end, and c) the port opens after the shock has reached 

the end. The constant pressure and velocity states that prevail in the 
passage just after the shock has reached the left end (time ’section* line 
in wave diagram), are perfectly realized in Fig. (I2a), while the contact 
surface is at the location shown by the sharp discontinuities in the density 
and entropy plots. Should the inlet port be opened earlier, e.g., at the time 
level shown by in the wave diagram, what happens is as follows: the 

pressure in the passage is still at the pre-compressed level and this comes 
into contact with the pressure level in the port which is considerably higher, 
resulting in' a shock propagating into the passage, colliding with the left 
moving shock and raising the overall pressure level to '^ 3.0 as shown in 
Fig. (I2b). However, as soon as the left moving shock reaches the end, it now 
encounters an open boundary with conditions that do not match those behind the 
shock, resulting in a rarefaction fan being generated, which propagates to the 
right. This expansion fan, travelling at sonic velocity relative to the gas 
into which it is propagating, soon overtakes the right moving shock which is 
travelling at a subsonic velocity relative to the same gas. This interaction 
results in an attenuation of both the rarefaction as well as the shock wave. 
Note that the overall pressure and velocity levels behind the rarefaction are 
about the same as for case a), i.e., the effects of the mismatch are not very 
significant at the outlet port. However, should the right moving pressure 
perturbations of case b) not attenuate each other significantly before they 
reach the right hand end, the consequences could be severe for the overall 
wave diagram, since this will lead to further (unwanted) wave reflections. 
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Fig. (12c) shows what occurs if the outlet port is opened too late, 
corresponding to time level on the wave diagram. Now the left 

travelling shock encounters a wall boundary condition on reaching the left 
end and reflects off as a shock, effectively doubling the pressure level 
behind it ( >3.5 in pressure plot of Fig. (12c)). l>/hen the outlet port opens, 
there is again a mismatch of conditions in the port and in the passage, with 
the pressure level in the passage being considerably higher than that 
prescribed for the outlet port. A rarefaction wave is generated which 
propagates to the right and overtakes the reflected shock. The same criterion 
holds for this case too, i.e., the ensuing attenuation of these pressure 
pulses should occur before they reach the right hand end, preferably even 
before they reach the interface still propagating towards the left at the 
flow velocity. 

The considerations above give a preview of the nature of decisions 
required in the successful design of a wave rotor device. It is clear that 
quite a few iterations are involved in the process of designing a viable wave 
diagram for a particular application, and each iteration entails calcul<iting 
two or more complete cycles to ensure ’closure’ or repeatability of the cycle. 
A fast solver like the RCM allows reaching an idealized ’base’ design quickly 
and inexpensively. 

Appendix A is a listing of the program in its present development 
stage. As mentioned earlier, no attempt has been made to optimize tlie 
program, either for storage requirements or for execution. 

Appendix B gives a description of the structure of the program, a 
listing of the important variables, the subroutines and the function 
subprograms. A step by step guide is also included to set up and run the 



program. 
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3. DISCUSSION AND RECOt^-lENDATIONS 



3.1. Discussion 

For meeting the criteria listed in the Introduction, in one dimension, 
Glimm*s method or the RCM appears to be superior to any difference method. 

For wave rotor type applications, where discontinuities need to be computed 
with sharpness, the 'infinite’ resolution of such discontinuities inherent in 
the RCM make it a natural choice to carry out ideal flow calculations for 
preliminary design purposes* Boundary conditions can be implemented quite 
easily and do not require information from points outside the domain of 
dependance as is the case in some finite difference schemes. The van der 
Corput sampling technique results in the best possible representation of the 
wave propagation, which is essential for the correct representation of 
continuous waves, particularly those produced by noalinear interactions. 

The method, however, is not recommended to solve for flows with real 
effects such as friction, heat transfer and area change, or to be extended to 
multi-dimensional flows. Although considerable research is being done to 
rigorously extend the method to such flows, with some degree of success (see 
Refs. 1, 4, 13), the present state of development is not mature enough to 
ensure a useful practical code as the outcome. 

3.2. Recoromendat ions 

Many options are available for one wishing to develop either a i-D code 
with real effects and/or a multi-dimensional code for wave rotor type 
applications. The author prefers to recommend numerical formulations which 
are dependent on the solution of Riemann problems, such as the Godunov method; 
the motivating reason for this preference is that a Riemann problem 
constitutes the solution of a discontinuity in the flow in terms of other 
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discontinuities (if any are, indeed, present), and the scheme is thus 
intrinsically suited for solving such flows; on the other hand, the other 
schemes, in general, require to be made aware of discontinuities in the flow 
through some external device, and then treat them through other artificial 
devices . 

A second-order, quasi one-dimensional (variable cross-sectional area) 
scheme has recently been developed by Ben-Artzi and Falcovitz (Ref, 14), The 
method is based on the exact solution of ’generalized Riemann propblems ’ , and 
has demonstrated very good results; it’s least accurate approximation is 
equivalent to Godunov’s first order method (Ref, 9), The resolution of shocks 
and other disconit inuities and singularities of the flow field is also high. 
Extension to more than one dimension appears to be straightforward through the 
use of operator splitting techniques, but has as yet not been tried 
extensively. 



25 



LIST OF REFERE.NCES 



1. Colella, P. , '*Gliram*s Method for Gasdynamics , SIAM Journal for 
Scientific Statistical Computations , Vol. 3, No, 1, March 1982, 

2. Glimm, J, , "Solutions in the Large for Nonlinear, Hyperbolic Systems 
of Equations," Communications in Pure and Applied Mathematics, No. 18, 
pp. 167-715, 19"6T: 

3. Chorin, A. J., "Random Choice Solution of Hyperbolic Systems," Journal of 
Computational Physics , No. 22, pp. 517-533, 1976. 

4. Chorin, A. J., "Random Choice Methods with Applications to Reacting Gas 
Flow," Journal of Computational Physics , No. 25, pp. 253-272, 1977. 

5. Godunov, S. K. , "A Finite Difference Method for the Numerical Computation 
and Discontinuous Solutions of the Equations of Fluid Dynamics," Mat. 
Sbornik , 47, pp. 357-393, 1959. 

6. Courant, R. and Friedrichs, K. 0., Supersonic Flow and Shock Waves , 
Interscience Publishers Inc., New York, 1948 

7. Hammersley, J. M. and Handscorab, D. C. , Monte Carlo Methods , John Wiley 
and Sons, Inc.-, New York, pp. 31-36, 1975. 

8. Rudinger , G. , Wave Diagrams for Nonstationary Flow in Ducts , Van 
Nostrand, 1955. 

9. Sod, G. A., "A Survey of Several Finite Difference Metliods for Systems of 
Non-linear Hyperbolic Conservation Laws," Journal of Computational 
Physics , No. 27, pp. 1-31, 1978. 

10. Shreeve, R. , Mathur , A., Eidelman, S. , and Erwin, J., Wave Rotor 
Technology Status and Research Progress Report , NPS 67-82-1) 14 PR, 
Turbopropulsion Laboratory, Naval Postgraduate School, Monterey, 
California, Novermber 1982. 

11. Klapproth, J. F. , Supercharged Turbowave Engines , General Electric 
R62FPD171, General Electric Flight Propulsion Division, Everdale, Ohio, 
May 1962. 

12. Taussig, R. , "Wave Rotor Turbofan Engines for Aircraft," Mechanical 
Engineering , Vol. 106, No. 11, pp. 60-68, November 1984. 

13. Gottlieb, J. J. and Igra, 0., "Interaction of Rarefaction Waves with Area 
Reductions in Ducts," Journal of Fluid Mechanics, Vol. 137, pp, 285-305, 
1983. 

14. Ben Artzi, M. and Falcovltz, J. , "A Second-Order Godunov Type Sclieme for 
Compressible Fluid Dynamics," Journal of Computational Physic s, No. 55, 
pp. 1-32, 1984. 



26 




FlG.l: THE RIC-IANN PRDBLEI^ IN GAS DITWUCS 
SPECLAL CASE OF SHOCK-TUBE FiaV 
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General Wave Sc rue Cure Resulting From Che Solution 
of a Riemann Problem 




Fig. 4 ; Ido;il Wave Diagram for Pre^buro Exchanger 
(Spectra Technology). 

I.P,IP,HP; Low, Intermediate, Mlgii Pressure 
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Test Case for 1-D Random Choice Method 
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Wave Diagram Computed by 1-D Random Choice 
Method. S — Shock; RS — Reflected Shock; 

R — ^refaction Pan; RR — Reflected Rarefaction; 
I — Interface; 

figure 6^ 
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Fig. 7 : Ideal Wave Diagram for General Electric 
Wave Engine 
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Fig. 8 : Gas Exhaust and Fresh Air Induction 
Process in G.E. Wave Engine 
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Fig. 10 a : Distribution of Flow Parameters in Fig. 10 b ; Distribution of Flow Parameters in Rotor 

Rotor Passage Just as Inlet Port Passage Just as Exhaust Port is Closed. 

Opens. N= 793. N= 998. 
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Fig. 10 c : Distribution of Flow Parameters in Fig. 11 : Distribution of Flow Parameters in Rotor 

Rotor Passage Just as Inlet Port Passage l^[hen Inlet Port Closes Late, i.e 

Closes. N= 1617. After Arrival Of Shock. N= 1700. 




Fig. 12 a : Distribution of Flow Parameters in 
Rotor Passage when the H.P. Air 
Outlet Port Opens on Time. 
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Fig. 12 b,c : Distribution of Flow Parameters in Rotor Passage 
for Two Conditions of Mismatch 
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Listing of Program RCM 
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PROGRAM ROM WITH VAN DER CORPUT SAMPLING AND SINGLE TIME STEP 
INTEGER QPRINT, QSTOP , SWL, SWR 
DIMENSION XX(6hYY(6) 

DIMENSION XARRAY(IOO) 

DIMENSION WNORM(12) , IDIGT(12) 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 

COMMON/ SUBS /P,R,U,A,S,X 

COMMON / GLIMMl / PGLIM , RGLIM , UGLIM , PL , RL , UL , PR , RR , UR , AL , AR , GL , GR , EPS 

COMMON / GLIMM2 / DT , DX , XI 

COMMON/ FUNl / G , PA , RA , UA , RB , RMU 

COMMON/ SAMPLE/ WNORM , IDIGT 

COMMON XARRAY,N1 

CALL COMPRS 

CALL BLOWUP (0.5) 

CALL PAGE (11. 0,8. 5) 

CALL HWSCAL( ’ SCREEN ' ) 

DATA K, SWL, SWR/ 500, 1,3/ 

DATA N , CFLNUM , TTOTAL/ 0,0.60,0.0/ 

DATA PSEXIT,PSINL,PSINR,RINL,RINR/116954. ,3770000. ,3819952.50,6.8 
*6.800/ 

DATA PSOUTl,PSOUT2,PSOUT3/ 38 19952. 5, 243 1800.0, 1530007 .5/ 

DATA PTOTIN,RTOTIN/ 1656663 .8,7.486/ 

DATA PREF,RREF,XREF/1656663.8,7.486,0. 1800/ 

G=1.4 

GL=1.4 

GR=1.4 

EPS=l.E-06 

QSTOP=20 

N1 = 0 

JC0UNT=0 

KCOUNT=0 

UEXMAX=0. 

DX=0.01 

AREF = SQRT ( PREF / RREF ) 

TIMREF=XREF/AREF 
RMU=SQRT( (G-1. )/ (G+1. ) ) 

X(1)=-0.5*DX 

ZETA=WDP(1) 

XII=DX*(WDP(0)-0.5) 

DO 25 1=2,203 
X(I)=X(I-1)+0.5*DX 
25 CONTINUE 

DO 35 1=1,100 
XARRAY(I)=X(I*2+1) 

35 CONTINUE 

INITIAL DATA 

CALL INITl 

CALL INIT2L(PSEXIT) 

CALL INIT2R(PSEXIT, PREF, RREF) 

CALL INIT3L(PSINL,RINL) 

CALL INIT3R(PSINR,RINR) 

NONDIMENSIONALIZATION 
DO 30 1=1,203,2 
P(I)=P(I)/PREF 



RCM00030 

RCM00040 

RCM00050 

RCM00060 

RCM00070 

RCM00080 

RCM00090 

RCMOOlOO 

RCMOOllO 

RCM00120 

RCM00130 

RCM00140 

RCM00150 

RCM00160 

RCM00170 

RCM00180 

RCM00190 

RCM00200 

RCM00210 

RCM00220 

RCM00230 

RCM00240 

RCM00250 

RCM00260 

RCM00270 

RCM00280 

RCM00290 

RCM00300 

RCM00310 

RCM00320 

RCM00330 

RCM00340 

RCM00350 

RCM00360 

RCM00370 

RCM00380 

RCM00390 

RCM00400 

RCM00410 

RCM00420 

RCM00430 

RCM00440 

RCM00450 

RCM00460 

RCM00470 

RCM00480 

RCM00490 

RCM00500 

RCM00510 

RCM00520 

RCM00530 

RCM00540 

RCM00550 

RCM00560 
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R(I)=R(I)/RREF RCMOi 

U(I)=U(I)/AREF RCMCI 

AU)=A(I)/AREF RCMC: 

S U ) = ALOG ( P ( I ) / R ( I ) — G ) RCMC 

30 CONTINUE RCMCI 

CALL PLOTl(K) RCMC! 

DO 40 J=1,K RCMC! 

N=N+1 RCMC 

XII=DX*(WDP(0)-0.5) RCMCi 

QPRINT=N/50 RCMC,* 

DT=100. RCMCI 

DO 50 1=1,203,2 RCMCI 

DTT=CFLNUM*DX/ (2.*AMAX1(ABS(A(I)+U(I) ) , ABS ( A( I ) -U ( I ) ) ) ) RCMCj^ 

DT=AMIN1(DTT,DT) RCMCjl 

50 CONTINUE RCMCi] 

TTOTAL=TTOTAL+DT RCMO^I 

TIME=TTOTAL*TIMREF RCMCI 

XI = -XII RCMCI 

DO 60 1 = 1,201,2 RCMCI 

PL=P(I) RCMC.] 

RL=R(I) RCMCI 

UL=U(I) RCMCj] 

PR=P(I + 2) RCMCj] 

RR=R(I + 2) . RCMCI 

UR=U(I + 2) RCMCI 

XITEMP=XI RCMCI 

IFU-EQ.1) XI = ABS(XI) RCMCI 

IF((I.EQ.201).AND. (XI.GT.0.0)) XI=-XI RCMC, I 

CALL GLIMM(QSTOP,PSTAR,USTAR,ASTAR) RCMCI 

XI = XITEMP RCMCI 

P(I+1)=PGLIM RCMCi 

R(I + 1)=RGLIM RCMCI 

UU+1)=UGLIM RCMCi 

60 CONTINUE RCMCi 

DO 70 1=1,201,2 RCMCj! 

IF(XI.LT.O.) GOTO 80 RCMC' 

P(I+2)=P(I+1) RCMC 

R(I + 2) = R(I+1) RCMC’' 

U(I-^2)=U(I + 1) RCMCI 

A(I + 2) = SQRT(G*P(I + 2)/R(I + 2)) RCMC:' 

S(I + 2)=AL0G(P(I + 2)/R(I + 2)-'"G) RCMC 

GOTO 70 RCMC ' 

80 P(I) = P(I + 1) RCMO' 

R(I) = RU-^1) RCMC 

U(I)=U(I+1) RCMC 

A(I) = SQRT(G'*P(I)/R(I)) RCMC 

S ( I ) =ALOG (P ( I ) /R( I ) — G ) RCMC 

70 CONTINUE RCMC. 

C CALL GE(SWL,SWR,N,TTOTAL,TIME,UEXMAX,PTOTIN,PREF) RCMC. 

C CALL DETON (SWL,SWR,N,QPRINT,TTOTAL, TIME) RCMC. 

CALL SPCTRA(N , SWL , SWR , TIME ,UEXMAX , PSEXIT , PSOUTl , PS0UT2 , PS0UT3 , J RCMC. 

*COUNT , QPRINT , TTOTAL , KCOUNT ) RCMC. 

IF(SWL.EQ.l) CALL BCLl RCMCj 

IF(SWL.EQ.2) CALL BCL2 (PSEXIT , PREF ) RCMCI 
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IF(SWL.EQ.3) CALL BCL3 ( PSINL , RINL , PREF , RREF ) RCMOlllC 

IF(SWL.EQ.4) CALL BCL4 (PTOTIN , RTOTIN , PREF , RREF ) * RCM0112C 

IF(SWL.EQ.5) CALL BCL5 RCM0113C 

IF(SWR.EQ.l) CALL BCR 1 RCM0114C 

IF(SWR.EQ.2) CALL BCR2 (PSEXIT , PREF ) RCM0115C 

IF(SWR.EQ.3) CALL BCR3 ( PSINR , RINR , PREF , RREF ) RCM0116C 

IF(SWR.EQ.4) CALL BCR4 ( PTOTIN , RTOTIN , PREF , RREF ) RCM0117C 

IF(SWR.EQ.5) CALL BCR5 RCM0118C 

IF ( (N.EQ. (50*QPRINT) ) .AND. (N.GE.O) ) CALL PLOT2(N,K) RCM0119C 

40 CONTINUE RCM0120C 

CALL ENDPL(O) RCM0121C 

CALL DONEPL RCM0122C 

STOP RCM0123C 

END RCM0124C 

SUBROUTINE GLIMM(QSTOP , PSTAR ,USTAR , ASTAR) RCM0125C 

INTEGER Q,QSTOP RCM0126G 

REAL ML,MR,MLN,MRN ' RCM0127C 

C0MM0N/GLIMM1/PGLIM,RGLIM,UGLIM,PL,RL,UL,PR,RR,UR,AL,AR,GL,GR,EPS RCM0128C 
COMMON /GLIMM2/DT,DX, XI RCM0129C 

DATA Q,ML,MR/0,100. ,100./ RCM0130C 

PSTAR=0.5*(PL+PR) RCM0131C 

COEFL=SQRT(PL''-RL) RCM0132C 

COEFR=SQRTiPR-RR) RCM0133C 

ALPHA=1. RCM0134G 

BEGIN GODUNOV ITERATION RCM0135G 

30 Q=Q+1 RCM0136G 

IF(PSTAR.LT.EPS) PSTAR=EPS RCM0137G 

COMPUTE NEXT ITERATION FOR ML AND MR RCM0138G 

MLN=COEFL*PHI (PSTAR, PL) RCM0139C 

MRN=COEFR"PHI (PSTAR, PR) RCM0140G 

DIFML=ABS(MLN-ML) RCM0141G 

DIFMR=ABS(MRN-MR) RCM0142G 

ML=MLN RCM0143G 

MR=MRN RCM0144G 

COMPUTE NEW PSTAR RCM0145G 

PTIL=PSTAR RCM01460 

PSTAR= (UL-UR+PL/ML+PR/MR)/ (1. /ML+1. /MR) RCM0147C 

PSTAR=ALPHA*PSTAR+ ( 1 . - ALPHA ) -PTIL RCM0148G 

IF(Q.LE.QSTOP) GOTO 10 RCM0149G 

IF(ABS(PSTAR-PTIL) .LT.EPS) GOTO 20 RCM0150C 

COMPUTE NEW ALPHA RCM0151G 

ALPHA=0. 5--ALPHA RCM0152G 

Q=0 RCM0153G 

IF((1. -ALPHA) .LT.EPS) GOTO 20 RCM0154C 

10 IF (DIFML.GE.EPS) GOTO 30 RCM0155C 

IF (DIFMR.GE.EPS) GOTO 30 RCM0156C 

END OF GODUNOV ITERATION; COMPUTE USTAR RCM0157C 

20 USTAR=(PL-PR+ML"UL+MR*UR)/(ML+MR) RCM0158G 

BEGIN SAMPLING PROCEDURE RCM0159C 

IF ( XI. LT. USTAR" DT) GO TO 40 RCM0160G 

RIGHT SIDE; SELECT CASE OF SHOCK OR EXPANSION RCM0161G 

IF (PSTAR.lt. PR) GO TO 50 RCM0162G 

RIGHT WAVE IS A SHOCK WAVE RCM0163C 

WR=UR+MR/RR RCM0164C 
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IF (XI.LT.WR“DT) GO TO 60 

C RIGHT OF RIGHT SHOCK CASE 

RGLIM=RR 
PGLIM=PR 
UGLIM=UR 
RETURN 

C LEFT OF RIGHT SHOCK CASE 

60 RGLIM=-MR/ (USTAR-WR) 

PGLIM=PSTAR 

UGLIM=USTAR 

RETURN 

C RIGHT WAVE IS A RAREFACTION WAVE 

50 CONST=PR/RR—GR 

RSTAR= ( PSTAR/ CONST )**( 1. /GR) 

ASTAR= SQRT ( GR*P STAR/ RSTAR ) 

AR = SQRT ( GR " PR / RR ) 

IF (XI.GE. (USTAR+ASTAR)*DT) GO TO 70 

C LEFT OF RIGHT FAN CASE 

RGLIM=RSTAR 
UGLIM=USTAR 
PGLIM=PSTAR 
RETURN 

C SELECT RIGHT OF FAN OR IN FAN 

70 IF (XI .GE. (UR+AR)*DT) GO TO 80 

C IN RIGHT FAN CASE 

UGLIM=2. / (GR+1. )*(XI/DT-AR+ (GR-1. )/2.*UR) 

RGLIM=( ( (AR+ (GR-1. ) / 2 . * (UGLIM-UR) )'"-2. )/ (GR-CONST) ) 

PGLIM= CONST-'RGLIM— GR 

RETURN 

C RIGHT OF RIGHT FAN CASE 

80 RGLIM=RR 
PGLIM=PR 
UGLIM=UR 
RETURN 

C LEFT SIDE; SELECT CASE OF SHOCK OR RAREFACTION 

40 IF (PSTAR.LT.pl) GO TO 90 

C LEFT WAVE IS A SHOCK WAVE 

WL=UL-ML/RL 

IF (XI.GE.WL*DT) GO TO 100 

C LEFT OF LEFT SHOCK CASE 

RGLIM=RL 
PGLIM=PL 
UGLIM=UL 
RETURN 

C RIGHT OF LEFT SHOCK CASE 

100 RGLIM=ML/ (USTAR-WL) 

PGLIM=PSTAR 

UGLIM=USTAR 

RETURN 

C LEFT WAVE IS A RAREFACTION WAVE 

90 CONST = PL/RL— GL 

RSTAR= (PSTAR/ CONST )**(!. /GL) 

ASTAR= SQRT ( GL*PSTAR / RSTAR ) 

AL= SQRT ( GL '• PL/ RL ) 
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IF (XI.LT. (USTAR-ASTAR)'"DT) GO TO 110 RCM02190 

RIGHT OF LEFT FAN CASE RCM0220C 

RGLIM=RSTAR RCM0221C 

PGLIM=PSTAR RCM02220 

UGLIM=USTAR RCM02230 

RETURN RCMO22A0 

SELECT LEFT OF FAN OR IN FAN CASE RCM02250 

110 IF (XI.LT. (UL-AL)*DT) GO TO 120 RCM02260 

IN LEFT FAN CASE RCM02270 

UGLIM=2. / (GL+1. )*(AL+ (GL-1. ) /2 . *UL+XI/DT) RCM02280 

RGLIM= ( ( (AL+(GL-1. ) / 2 . - (UL-UGLIM) )-'--'2. )/ (GL--CONST) )"*(1. / (GL-1. )) RCM02290 
PGLIM=CONST*RGLIM““GL RCM02300 

RETURN RCM02310 

LEFT OF LEFT FAN CASE RCM02320 

120 RGLIM=RL RCM02330 

PGLIM=PL RCM02340 

UGLIM=UL RCM02350 

RETURN RCM02360 

END RCM02370 

FUNCTION PHI(Y,Z) RCM02380 

REAL RMU RCM02390 

COMMON/ FUN 1 / G , PA , RA , UA , RB , RMU RCMO 24 0 0 

EPS=l.E-06- RCM02410 

PARAM=Y/Z RCM02420 

IF (ABS(1. -PARAM) .GE.EPS) GO TO 10 RCM0243C 

PHI=SQRT(G) RCM02440 

RETURN RCM02450 

10 IF (PARAM. GE. 1. ) GO TO 20 RCM02460 

PHI= (G- 1 . ) / 2 . * ( 1 . -PARAM) / ( SQRT(G)*( 1 . - PARAM— (( G- 1 . ) / ( 2 . - G ) ) ) ) RCM02470 

RETURN RCM02480 

20 PHI=SQRT((G+1. )/2."PARAM+(G-l. )/2. ) RCM02490 

RETURN RCM02500 

END RCM02510 

FUNCTION PHIl(PB) RCM02520 

REAL RMU RCM02530 

COMMON/FUNl/G,PA,RA,UA,RB ,RMU RCM02540 

PHI1= (PB-PA)"'SQRT( (1. -RMU*~2. )/ (RA« (PB + RMU"*2 . *PA) ) ) RCM02550 

RETURN RCM02560 

END RCM02570 

FUNCTION PSI(PB) RCM02580 

REAL RMU RCM02590 

COMMON / FUN 1 / G , PA , RA , UA , RB , RMU RCMO 2600 

PSI = SQRT(1. -RMU-~4. ) /RMU — 2 . / SQRT (RA)*PA** ( 1 . / ( 2 . “G ) )* (PB-* -• ( (G- 1 . )RCM02610 
* / ( 2 . -G ) ) - PA-"' ( ( G - 1 . ) / ( 2 . -G ) ) ) RCM02 6 20 

RETURN RCM02630 

END RCM02640 

SUBROUTINE INITl RCM02650 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) RCM02660 

COMMON/ FUNl/G, PA, RA,UA,RB, RMU RCM02670 

COMMON/SUBS/P, R,U, A, S,X RCM02680 

DO 10 1=1, 9, 2 RCM02690 

P(I)=810600.00 RCM02700 

RU)=0.7132 RCM02710 

U(I)=644.4 RCM02720 
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A(I) = SQRT(G*P(I)/R(I) ) 

10 CONTINUE 

DO 20 1=11,203,2 
P(I)=101325.0 
R(I)=1.22 
U ( I ) = 0 . 0 

A(I)=SQRT(G*P(I)/R(I) ) 

20 CONTINUE 
RETURN 
END 

SUBROUTINE INIT2R ( PSEXIT , PREF , RREF ) 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) , S (203 ) , X (203 ) 

COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 

COMMON/ SUBS/ P,R,U, A, S,X 

DO 10 1=3,201,2 

P(I)=PREF 

R(I)=RREF 

U(I)=0.0 

A(I )=SQRT(P(I)*G/R(I) ) 

10 CONTINUE 
P(1)=P(3) 

R(1)=R(3) 

U(l) = -U(3) . 

A(1) = SQRT(G''P(1)/R(1) ) 

P(203)=PSEXIT 

R(203)=R(201) 

PA=P(201) 

RA=R(201) 

UA=U(201) 

PB=P(203) 

RB=R(203) 

IF(PA.GT.PB) GO TO 20 
U(203)=UA-PHI1(PB) 

GO TO 30 

20 U(203)=UA-PSI(PB) 

30 A(203) = SQRT(G*P(203)/R(203) ) 

RETURN 

END 

SUBROUTINE INIT2L(PSEXIT) 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) , S (203 ) ,X(203 ) 

COMMON / FUNl / G , PA , RA , UA , RB , RMU 

COMMON/ SUBS /P,R,U, A, S,X 

DO 10 1=3,201,2 

P(I)=285080.0 

R(I)=0.897 

U(I)=0.0 

AU) = SQRT(G*P(I)/R(I) ) 

10 CONTINUE 
RETURN 
END 

SUBROUTINE INIT3L(PSINL,RINL) 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 
COMMON/ FUNl / G , PA , RA , UA , RB , RMU 
COMMON/ SUBS / P , R , U , A , S , X 



DO 10 1=3,201,2 
P(I)=2390000.0 
R(I)=9.787 
U(I)=0.0 

A(I)=SQRT(G*P(I)/R(I) ) 

10 CONTINUE 
P(1)=PSINL 
R(1)=RINL 
PA=P(3 ) 

RA=R(3) 

UA=U(3) 

PB=P(1) 

U(1)=UA+PHI1(PB) 

A(1)=SQRT(G"P(1)/R(1) ) 

P(203)=P(201) 

R(203)=R(201) 

U(203)=-U(201) 

A(203) = SQRT(G-'P(203)/R(203) ) 

RETURN 

END 

SUBROUTINE INIT3R (PSINR , RINR) 

DIMENSION P(203) ,R(203) ,U(203) ,A(203 ) , S (203 ) ,X(203 ) 

COMMON/SUBS /P,R,U, A, S ,X 

COMMON/ FUNl/ G , PA , RA , UA , RB , RMU 

DO 10 1=3,201,2 

P(I)=2421667.5 

R(I)=9.787 

U(I)=0.0 

A(I) = SQRT(G“P(I)/R(I) ) 

10 CONTINUE 

P(203)=PSINR 

PA=P(201) 

RA=R(201) 

UA=U(201) 

PB=P(203) 

U(203)=UA-PHI1(PB) 

R(203)=RINR 

A(203)=SQRT(G*P(203)/R(203) ) 

P(1)=P(3) 

R(1)=R(3) 

U(l)=-U(3) 

A(1)=SQRT(G*P(1)/R(1) ) 

RETURN 

END 

SUBROUTINE BCLl 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 
COMMON/ SUBS /P,R,U, A, S,X 
COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 
P(1)=P(3) 

R(1)=R(3) 

U(l)=-U(3) 

A(1)=SQRT(G“P(1)/R(1) ) 

RETURN 

END 



RCM03270 

RCM03280 

RCM03290 

RCM03300 

RCM03310 

RCM03320 

RCM03330 

RCM03340 

RCM03350 

RCM03360 

RCM03370 

RCM03380 

RCM03390 

RCM03400 

RCM03410 

RCM03420 

RCM03430 

RCM03440 

RCM03450 

RCM03460 

RCM03470 

RCM03480 

RCM03490 

RCM03500 

RCM03510 

RCM03520 

RCM03530 

RCM03540 

RCM03550 

RCM03560 

RCM03570 

RCM03580 

RCM03590 

RCM03600 

RCM03610 

RCM03620 

RCM03630 

RCM03640 

RCM03650 

RCM03660 

RCM03670 

RCM03680 

RCM03690 

RCM03700 

RCM03710 

RCM03720 

RCM03730 

RCM03740 

RCM03750 

RCM03760 

RCM03770 

RCM03780 

RCM03790 

RCM03800 
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SUBROUTINE BCRl RCMO 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) RCMOi 

COMMON/SUBS/P, R,U, A, S,X RCMOj 

COMMON/FUNl/G,PA,RA,UA,RB,RMU RCMO) 

P(203)=P(201) RCMO, 

R(203)=R(201) RCMO 

U(203)=-U(201) RCMOi 

A(203)=SQRT(G*P(203)/R(203)) RCMOi 

RETURN RCMOi 

END RCMOi 

SUBROUTINE BCL2 (PSEXIT , PREF ) RCMO; 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) RCMOi 

COMMON/SUBS/P, R,U, A, S,X RCMO; 

COMMON/FUNl/G,PA,RA,UA,RB,RMU RCMO 

P(1) = PSEXIT/PREF RCMO' 

R(1)=R(3) RCMOI 

U(1)=U(3)- RCMOI 

A(1)=SQRT(G*P(1)/R(1)) RCMOj 

RETURN RCMO I 

END RCMO| 

SUBROUTINE BCR2 ( PSEXIT , PREF ) RCMO! 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) RCMOI 

COMMON/SUBS/P, R,U, A, S,X RCMO.'l 

COMMON/FUn1/G,PA,RA,UA,RB,RMU - RCMO.) 

P(203) = PSEXIT/PREF RCMO:) 

R(203)=R(201) RCMOa 

PA=P(20l) RCMOli 

RA=R(201) RCMO 13 

UA=U(201) RCMOI 3 

PB=P(203) RCMOI 

RB=R(203) RCMOI 

IF(PA.GT.PB) GO TO 10 RCMOI 

U(203)=UA-PHI1(PB) RCMOI 

GO TO 20 RCMOI 

10 U(203)=UA-PSI(PB) RCMOa 

20 A(203)=SQRT(G*P(203)/R(203) ) RCMOI 

RETURN RCMOI 

END RCMOP 

SUBROUTINE BCL3 (PSINL , RINL , PREF , RREF ) RCMOI 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) RCMO| 

COMMON/SUBS/P, R,U, A, S,X RCMO| 

C0MM0N/FUN1/G,PA,RA,UA,RB,RMU RCMOi 

P(1) = PSINL/PREF RCMO'j 

R(1)=RINL/RREF RCMOjj 

PA=P(3) RCMOI 

RA=R(3) RCMO; 

UA=U(3) RCMO| 

PB = P(1) RCMOi 

U(1)=UA+PHI1(PB) RCMO; 

A(1)=SQRT(G“P(1)/R(1)) RCMO 

RETURN RCMO 

END RCMO 

SUBROUTINE BCR3 (PSINR , RINR , PREF , RREF ) RCMO 

DIMENSION P(203),R(203) ,U(203),A(203),S(203),X(203) RCMO 

I 
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COMMON/ SUBS /P,R,U, A, S,X 


RCM04350 




COMMON/ FUNl / G , PA , RA , UA , RB , RMU 


RCM04360 




P(203)=PSINR/PREF 


RCM04370 




R(203)=RINR/RREF 


RCM04380 




PA=P(201) 


RCM04390 




RA=R(201) 


RCM04400 




UA=U(201) 


RCM04410 




PB=P(203) 


RCM04420 




U(203)=UA-PHI1(PB) 


RCM04430 




A(203)=SQRT(G*P(203)/R(203)) 


RCM04440 




RETURN 


RCM04450 




END 


RCM04460 




SUBROUTINE BCL4 (PTOTIN , RTOTIN , PREF , RREF ) 


RCM04470 




INTEGER QOUT 


RCM04480 




DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 


RCM04490 




DIMENSION XARRAYUOO) 


RCM04500 




COMMON/ SUBS /P,R,U, A, S,X 


RCM04510 




COMMON/ FUNl/ G , PA , RA , UA , RB , RMU 


RCM04520 




COMMON XARRAY,N1 


RCM04530 




N1=N1+1 


RCM04540 




QOUT=Nl/5 


RCM04550 




PTOT=PTOTIN/PREF 


RCM04560 




RTOT=RTOTIN-/RREF 


RCM04570 




ATOT=SQRT(G“PTOT/RTOT) 


RCM04580 




STOT=ALOG (PTOT/RTOT—G ) 


RCM04590 




U(1)=U(3) 


RCM04600 




A ( 1 ) = SQRT ( ATOT- 2 . - ( G- 1 . ) / 2 . * ABS (U ( 1 ) ) —2 . ) 


RCM04610 




AMACH=U(1)/A(1) 


RCM04620 




IF (AMACH.lt. 0.0) GO TO 60 


RCM04630 




P(l)=PTOT/(l.+(G-l. ) / 2 . *AMACH- *2 . )**(G/(G-1. )) 


RCM04640 




R(l)=RTOT/ (1.+ (G-1. )/2. "AMACH**2. )*-(l. / (G-1. ) ) 


RCM04650 




S ( 1) =ALOG ( P ( 1 ) /R( 1 )-" "G ) 


RCM04660 




GO TO 50 


RCM04670 


60 


P(l) = PTOT/ (1. + (G-1. )/2. "ABS(AMACH)— 2. )*-(G/ (G-1. ) ) 


RCM04680 




R(l) = RTOT/ (1. + (G-1. )/2. "'ABS(AMACH)*''2. ) '”• ( 1 . / ( G- 1 . )) 


RCM04690 




S ( 1 ) = ALO G ( P ( 1 ) / R ( 1 ) * " G ) 


RCM04700 


50 


RETURN 


RCM04710 




END 


RCM04720 




SUBROUTINE BCR4 ( PTOTIN , RTOTIN , PREF , RREF ) 


RCM04730 




INTEGER QOUT 


RCM04740 




DIMENSION P(203) ,R(203) ,U(203) ,A(203) , S ( 203 ) , X ( 203 ) 


RCM04750 




DIMENSION XARRAYUOO) 


RCM04760 




COMMON/ SUBS /P,R,U, A, S,X 


RCM04770 




COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 


RCM04780 




COMMON XARRAY,N1 


RCM04790 




N1=N1+1 


RCM04800 




QOUT=Nl/25 


RCM04810 




PTOT=PTOTIN/PREF 


RCM04820 




RTOT=RTOTIN/RREF 


RCM04830 




ATOT= SQRT ( G --'PTOT/ RTOT ) 


RCM04840 




STOT=ALOG(PTOT/RTOT**G) 


RCM04850 




U(203)=U(201) 


RCM04860 




A ( 20 3 ) = SQRT ( ATOT** 2 . - ( G - 1 . ) / 2 . *ABS (U ( 20 3 ) ) **2 . ) 


RCM04870 




AMACH=U(203)/A(203) 

49 
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IF (AMACH.lt. 0.0) GO TO 60 


RCM04 


P(203)=PTOT/(1. + (G-1. )/2. ~AMACH*“2. ) — (G/(G-1. )) 


RCMO^ 


R(203)=RTOT/(1.+(G-1. ) /2 . *AMACH*"2 . )— (l./(G-l. )) 


RCMO-!: 


SU03)=ALOG(P(203)/R(203)**G) 


RCMO^i 


GO TO 50 


RCMO^I 


60 P(203)=PTOT/(1.+(G-1. ) / 2 . -ABS (AMACH)**2 . ) — (G/ (G- 1 . )) 


RCMO^i 


R(203)=RTOT/ (1. + (G-1. ) /2 . -ABS (AMACHr~*2 . )**(!. /(G-1. )) 


RCMO^i 


S(203) = ALOG(P(203)/R(203)*’"G) 


RCMO^ 


50 RETURN 


RCMO^ 


END 


RCMO^ 


SUBROUTINE BCL5 


RCMO^! 


DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 


RCMO: 


COMMON/ SUBS /P,R,U, A, S,X 


RCMOI 


P(1)=P(3) 


RCMO.'' 


R(1)=R(3) 


RCM0I‘^ 


U(1)=U(3) 


RCMO!=; 


A(1)=A(3) 


RCMOi’l 


RETURN 


RCMO:*'l 


END 


RCMOlli 


SUBROUTINE BCR5 


RCMO^fl 


DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 


RCMOIjJ 


COMMON/ SUBS/P, R,U, A, S,X 


RCMOLI 


P(203)=P(201) 


RCMO![l 


R(203)=R(201) 


RCMO:'.! 


U(203)=U(201) 


RCMOij! 


A(203)=A(201) 


RCMO:j.l 


RETURN 


RCMO:.' 


END 


RCMO; ( 


FUNCTION WDP(II) 


RCMO!;j 


DIMENSION WNORM(12) ,IDIGT(12) 


RCMOy 


COMMON/ SAMPLE/ WNORM , IDIGT 


RCMO:.l 


IF (lI.EQ.O) GO TO 10 


RCMOlfl 


Ll=2 


RCMOJij 


L2 = l 


RCMO!,;; 


DO 20 JJ=1,12 


RCMO!^i 


IDIGT(JJ)=0 


RCMOiP 


V7NORM(JJ) = l. /FL0AT(L1 — JJ) 


RCMOli'i 


20 CONTINUE 


RCMO 21 


WDP=0. 


RCMO.f 


RETURN 


RCMO;J| 


10 DO 40 JJ=1,12 


RCMO;S 


Ll = 2 


RCMOU 


L2=l 


RCMO.® 


KJO=IDIGT(JJ) 


RCMO. a 


KJN=MOD( (KJO+1) ,L1) 


RCMO.il 


IDIGT (JJ)=KJN 


RCMO. a 


IF (KJO.LT.KJN) GO TO 50 


RCMO.i) 


40 CONTINUE 


RCMO.'; 


50 SUM=0. 


RCMO.j: 


DO 60 JJ=1,12 


RCMO.; 


KNEW=MOD ( IDIGT (JJ )*L2 , LI ) 


RCMO.i] 


SUM= SUM+ FLOAT (KNEW ) "WNORM ( J J ) 


RCMO. 5 


60 CONTINUE 


RCMO.li 


WDP=SUM 


RCMO.; 
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I 



RETURN 

END 

SUBROUTINE PLOTl(K) 

DIMENSION X0RG(4) ,YORG(A) ,YMAX(4) ,YMIN(4) 
DATA XORG/0.5,4.75,0.5,4.75/ 

DATA YORG/ 0.5, 0.5, 4. 75, 4. 75/ 

DATA YMAX/3.50,3.0,0.5,2.0/ 

DATA YMIN/ 0.5, 0.0, -0.5, -2.0/ 

DO 10 1=1,4 



CALL 

CALL 

CALL 

CALL 

CALL 



GRAF(0. , 
ENDGR(O) 
10 CONTINUE 

CALL PHYSOR(8 
AREA2D(2 
FRAME 
GRAF(0. , 
ENDGR(O) 



PHYSOR ( XORG ( I ) , YORG ( I ) ) 
AREA2D(3.5,3.5) 

FRAME 

SCALE’ ,1.0,YMIN(I) 



SCALE' ,YMAX(I) ) 



5,0.5) 

25,7.75) 



SCALE’ , 1. ,0, ’SCALE’ ,K) 



CALL 
CALL 
CALL 
CALL 
RETURN 
END 

SUBROUTINE. PLOT2 (N ,K) 

DIMENSION XORG(4) ,YORG(4) ,YMAX(4) ,YMIN(4) ,KNT(4 ) , IYNAM( 10 ) 
DIMENSION PARRAY(IOO) ,RARRAY(100) ,UARRAY(100) ,SARRAY(100) ,XARRAY( 
* 00 ) 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) , S (203 ) ,X(203 ) 

COMMON/ SUBS /P,R,U, A, S,X 
COMMON XARRAY 

DATA XORG/0.5,4.75,0.5,4.75/ 

DATA YORG/0.5,0.5,4.75,4.75/ 

DATA YMAX/3.50,3.0,0.5,2.0/ 

DATA YMIN/ 0.5, 0.0, -0.5, -2.0/ 

DATA KNT/1,4,6 ,9/ 

DATA lYNAM/ ’ PRES ’ , ’ SURE ’ , ’ $ 

*, ’ENTR' , ’OPY$ ’ / 

DO 200 1=1,100 
PARRAY(I)=P(I*2+1) 

RARRAY(I)=R(I*2+1) 

UARRAY(I)=U (1*2+1) 

SARRAYU) = SU''2 + 1) 

200 CONTINUE 

DO 300 1=1,4 

CALL PHYSOR(XORG(I) ,YORG(I) ) 

AREA2D(3.5,3.5) 

XNAME ( ’ X ’ , 1 ) 

YNAME ( lYNAM (KNT ( I ) ) , 100 ) 

GRAF(0. , ’SCALE’ ,1.0,YMIN(I) 

EQ. 1) CALL SETCLR( ’YELLOW’ ) 



DENS ’ , ’ ITY$ ’ , ’ VELO ’ , ’ CITY ’ , ’ $ 



CALL 
CALL 
CALL 
CALL 
IF(I. 
IF(I.EQ.2) 
IF(I.EQ.3) 
IF(I.EQ.4) 
IF(N.EQ.K) 
IFU .EQ. 1) 



SCALE’ , YMAX(I) ) 



CALL SETCLR( ’CYAN’ ) 

CALL SETCLR( ’RED’ ) 

CALL SETCLR ( ’ MAGENTA ’ ) 

CALL SETCLR ( ’WHITE’ ) 

CALL CURVE (XARRAY , P ARRAY , 100 , 0 ) 



RCM05430 

RCM05440 

RCM05450 

RCM05460 

RCM05470 

RCM05480 

RCM05490 

RCM05500 

RCM05510 

RCM05520 

RCM05530 

RCM05540 

RCM05550 

RCM05560 

RCM05570 

RCM05580 

RCM05590 

RCM05600 

RCM05610 

RCM05620 

RCM05630 

RCM05640 

RCM05650 

RCM05660 

1RCM05670 

RCM05680 

RCM05690 

RCM05700 

RCM05710 

RCM05720 

RCM05730 

RCM05740 

RCM05750 

RCM05760 

’RCM05770 

RCM05780 

RCM05790 

RCM05800 

RCM05810 

RCM05820 

RCM05830 

RCM05840 

RCM05850 

RCM05860 

RCM05870 

RCM05880 

RCM05890 

RCM05900 

RCM05910 

RCM05920 

RCM05930 

RCM05940 

RCM05950 

RCM05960 
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I 



( XARRAY , RARRAY ,100,0) 
(XARRAY , UARRAY ,100,0) 
(XARRAY , SARRAY , 100 , 0 ) 



GO TO 10 
GO TO 30 
GO TO 50 
RETURN 



10 



20 



40 



50 



60 



) 



RCM(; 

RCM(j 

RCM(tj 

RCM([ 

RCM(^ 

RCM( 

RCM( 

RCM( 

RCM( 

RCM( 

RCM( 

RCM( 

RCM( 

RCM( 

RCM( 

RCM( 

RCM( 

RCM( 

RCM(i 

RCM(k 

RCM( 

RCMC 



RCM( 



IF(I.EQ.2) CALL CURVE 
IF(I.EQ.3) CALL CURVE 
IF(I.EQ.4) CALL CURVE 
CALL ENDGR(O) 

300 CONTINUE 
RETURN 
END 

SUBROUTINE GE ( SWL , SWR , N , TTOTAL , TIME , UEXMAX , PTOTIN , PREF ) 

INTEGER SWL, SWR 

DIMENSION P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 

COMMON/ SUBS /P,R,U,A,S,X 

C — CALCULATION STARTS AT EXHAUST PORT OPENING. SUBROUTINE STRUCTURED 
C**ACCORDINGLY . 

IF((SWL.EQ. 1) .AND. (SWR.EQ.2)) 

IF ( ( SWL . EQ . 4 ) . AND . ( SWR . EQ . 2 ^ 

IF ( ( SWL . EQ . 4 ) . AND . ( SWR . EQ . 1 ) ) 

IF ( ( SWL . EQ . 1 ) . AND . ( SWR . EQ . 1 ) ) 

PWALL=P(2) 

IF(PWALL.LE. (PTOTIN/ PREF) ) GO TO 20 
RETURN 
SWL = 4 

WRITE (6, 74) 

WRITE (6, 7 5J N, TTOTAL, TIME 
RETURN 

30 UEXIT=U(202) 

IF (UEXMAX. LT.UEXIT) UEXMAX=UEXIT 
IF (UEXIT.lt. UEXMAX/ 2. ) GO TO 40 
RETURN 
SWR=1 

WRITE (6, 76) 

WRITE (6, 75) N, TTOTAL, TIME 
RETURN 
P1SH0K=P(2) 

IF(P1SH0K.GT.PT0TIN/PREF) GO TO 60 
RETURN 
SWL=1 

WRITE (6, 77) 

WRITE (6,75) N , TTOTAL , TIME 

74 FORMAT ( 5 X, 'INLET PORT OPENS AT:') 

75 F0RMAT(5X,I4,5X,2F14. 7) 

76 FORMAT (5X, 'EXHAUST PORT CLOSES AT: 

77 FORMAT ( 5X, ' INLET PORT CLOSES AT:') 

RETURN 
END 

SUBROUTINE SPCTRA(N , SWL , SWR , TIME , UEXMAX , PSEXIT , PSOUTl , PSOUT2 , PSOUTRCM(j) 
*3 ,JCOUNT,QPRINT, TTOTAL, KCOUNT) RCM(» 

INTEGER SWL,SWR,QPRINT RCM(ji 

DIMENSION pU 03) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) RCM(i 

COMMON/SUBS/P, R,U, A, S ,X RCMC 

C*CALCULATION STARTS AT HP GAS IN PORT. JCOUNT IS NUMBERED ACCORDINGLY**RCM( 
IF((SWL.EQ. 1) .AND. (SWR.EQ.3) ) GO TO 10 RCMC 

IF((SWL.EQ.2) .AND. (SWR.EQ.3)) GO TO 20 RCM(|i 

IF((SWL.EQ.2) .AND. (SWR.EQ.l)) GO TO 30 RCMC 

IF( (SWL.EQ. 5) .AND. (SWR.EQ. 1) ) GO TO 40 RCMCi 



RCMC 

RCMC 

RCMC 

RCMC 

RCMC 

RCMCi, 

RCMCl 

RCMC! 

RCMC! 

RCMC 

RCMC 

RCMC 

RCMC) 

RCMC* 

RCMO 

RCMC. 

RCMC). 

RCM(. 

RCMCI) 

RCMCf 

RCMCi) 
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IF((SWL.EQ.2) .AND. (SWR.EQ.5)) GO TO 50 
IF((SWL.EQ.2).AND.(SWR.EQ.4^ GO TO 60 
IF((SWL.EQ. 1) .AND. (SWR.EQ.4h GO TO 70 

10 IF(U(3) .LT.0.0) GO TO 11 
RETURN 

11 JC0UNT=JC0UNT+1 
PSEXIT=PS0UT1 
SWL=2 

WRITE(6 , 12) 

WRITE(6 , 13)N,TIME,SWL,SWR, JCOUNT 

12 FORMAT ( 5X, 'H.P. AIR OUT PORT OPENS AT') 

13 F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

20 DO 26 1=5,199,2 

IF((R(I)-R(I+2) ) .GT.O. 1) GO TO 22 
GO TO 26 

22 XCNTCT=X(I) 

UCNTCT=U(I) 

TCNTCT=XCNTCT/ABS(UCNTCT) 

AHEAD=A(199) 

THEAD=1.0/AU99) 

IF(TCNTCT.LE.THEAD) GO TO 23 
RETURN 

23 JCOUNT=JCOUNT+l 
SWR=1 

WRITE(6 ,24) 

WRITE ( 6 , 2 5 )N , TIME , SWL , SWR , JCOUNT 

24 FORMAT (5X, 'H.P. GAS IN PORT CLOSES AT') 

25 F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

26 CONTINUE 

30 IF ( JCOUNT. EQ. 4) GO TO 80 
IF ( JCOUNT. EQ. 6) GO TO 90 
IF ( JCOUNT. EQ. 8) GO TO 100 

IF( (R(2)-R(4) ) .GT.O. 1) GO TO 31 
RETURN 

31 JCOUNT=JCOUNT+l 
SWL=5 

WRITE(6,32) 

WRITE ( 6 , 3 3 )N , TIME , SWL , SWR , JCOUNT 

32 FORMAT(5X, 'HP AIR OUT PORT CLOSES AND TUNING PORT LI OPENS AT') 

33 F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

40 IF(JCOUNT.EQ.7) GO TO 110 
IF(U(3) .GE.0.0) GO TO 41 
RETURN 

41 JCOUNT=JCOUNT+l 
PSEXIT=PSOUT2 
SWL=2 

WRITE(6,42) 

WRITE(6 ,43)N,TIME,SWL,SWR, JCOUNT 

42 FORMAT (5X, 'TUNING PORT LI CLOSES AND EXHAUST PORT El OPENS AT') 

43 F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 



RCM0651C 
RCM0652C 
RCM0653C 
RCM0654C 
RCM0655C 
RCM0656C 
RCM0657C 
RCM0658C 
RCM0659C 
RCM0660C 
RCM0661C 
RCM06620 
RCM0663C 
RCM06640 
RCM06650 
RCM06660 
RCM06670 
RCM06680 
RCM06690 
RCM06700 
RCM06710 
RCM06720 
RCM06730 
RCM0674C 
RCM0675C 
RCM06760 
RCM06770 
RCM06780 
RCM06790 
RCM06800 
RCM068 10 
RCM06820 
RCM06830 
RCM06840 
RCM06850 
RCM06860 
RCM06870 
RCM06880 
RCM06890 
RCM06900 
RCM06910 
RCM06920 
RCM06930 
RCM06940 
RCM06950 
RCM06960 
RCM06970 
RCM06980 
RCM06990 
RCM07000 
RCM07010 
RCM07020 
RCM07030 
RCM07040 



. ti 
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80 IF(ABS(U(201)) ,GT. .0001) GO TO 81 
RETURN 

81 JCOUNT=JCOUNT+1 
SWR=5 

WRITE(6 ,82) 

WRITE ( 6 , 8 3 ) N , TIME , SWL , SWR , JCOUNT 

82 FORMAT ( 5 X, 'TUNING PORT R1 OPENS AT') 

83 F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

50 IF(JCOUNT.EQ. 9) GO TO 120 
IF(ABS(U(201)-U(3)) .LE. 0.001) GO TO 51 
RETURN 

51 JCOUNT=JCOUNT+l 
SWR=1 

WRITE(6 ,52) 

WRITE ( 6 , 5 3 )N , TIME , SWL , SWR , JCOUNT 

52 FORMAT ( 5 X, 'TUNING PORT R1 CLOSES AT') 

53 F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

90 THEAD1=X(201)/(ABS(U(3))+A(3)) 

KCOUNT=KCOUNT+l 
IF(KCOUNT.EQ. 1) TTOTl=TTOTAL 
IF(TTOTAL.GE. (TTOTl+THEADl) ) GO TO 91 
RETURN 

91 JCOUNT=JCOUNT+l 
SWL=5 

WRITE(6 ,92) 

WRITE ( 6 , 9 3 ) N , TIME , SWL , SWR , JCOUNT 

92 FORMAT (5X, 'EXHAUST PORT El CLOSES AND TUNING PORT L2 OPENS AT') 

93 F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

110 IF(U(3) .GE. 0.0) GO TO 111 
RETURN 

111 JCOUNT=JCOUNT+l 
PSEXIT=PSOUT3 
SWL=2 

WRITE(6,112) 

WRITE ( 6 , 1 13 ) N , TIME , SWL , SWR , JCOUNT 

112 FORMAT (5X, 'TUNING PORT L2 CLOSES AND EXHAUST PORT E2 OPENS AT') 

113 F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

100 IF(ABS(U(201) ) .GT. .0001) GO TO 101 
RETURN 

101 JCOUNT=JCOUNT+l 
SWR=5 

WRITE (6, 102) 

WRITE ( 6 , 103 )N , TIME , SWL , SWR , JCOUNT 

102 FORMAT (5X, 'TUNING PORT R2 OPENS AT') 

103 F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

120 IF(ABS(U(201)-U(3) ) .LE.O .0001) GO TO 121 
RETURN 

121 JCOUNT=JCOUNT+l 
SWR=4 



RCM07(» 

RCM07(| 

RCM07(| 

RCM07(i 

RCM07(ll 

RCM07;|| 

RCM07;| 

RCM07;i 

RCM07; 

RCM07: 

RCM07; 

RCM07: 

RCM07; 

RCM07:| 

RCM07:'| 

RCM07;| 

RCM07:| 

RCM07:' 

RCM07:; 

RCM07:' 

RCM07:: 

RCM07: 

RCM07: 

RCM07.'i 

RCM07,": 

RCM07.| 

RCM07i) 

RCM07:ii 

RCM07.i 

RCM07!I 

RCM073 

RCM07i< 

RCM07^‘ 

RCM07: 

RCM07 ' 

RCM07', 

RCM07j 

RCM07j 

RCM07I 

RCM07 

RCM07 

RCM07| 

RCM07| 

RCM07! 

RCM07; 

RCM07- 

RCM07j 

RCM07: 

RCM07j 

RCM07: 

RCM07 

RCM07 

RCM07 

RCM07 
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WRITE(6, 122) 

WRITE(6 , 123 ) N, TIME, SWL,SWR,JCOUNT 

FORMAT ( 5 X, ’TUNING PORT R2 CLOSES AND L.P. AIR INLET OPENS AT’) 


RCM07590 

RCM07600 


122 


RCM07610 


123 


FORMAT ( 5X , 14 , 5X , F9 . 7 , 5X , 313 ) 
RETURN 


RCM07620 

RCM07630 


60 


IF( (R(4)-R(2) ) .GT.O . 1) GO TO 61 
RETURN 


RCM07640 

RCM07650 


61 


JCOUNT=JCOUNT+l 

SWL=1 

WRITE (6, 62) 

WRITE ( 6 , 6 3 ) N , TIME , SWL , SWR , JCOUNT 
FORMAT (5X, ’EXHAUST PORT E2 CLOSES AT’) 


RCM07660 

RCM07670 

RCM07680 

RCM07690 


62 


RCM07700 


63 


F0RMAT(5X,I4,5X,F9. 7 ,5X,3I3) 
RETURN 


RCM07710 

RCM07720 


70 


IF(U(201) .GE.0.0) GO TO 71 
RETURN 


RCM07730 

RCM07740 


71 


JC0UNT=0 

SWR=1 

WRITE (6 ,72) 

WRITE ( 6 , 7 3 ) N , TIME , SWL , SWR , JCOUNT 
FORMAT (5X, ’CYCLE COMPLETED.’) 


RCM07750 

RCM07760 

RCM07770 

RCM07780 


72 


RCM07790 


73 


FORMAT ( 5X , 14 , 5X , F9 . 7 , 5X , 313 ) 

RETURN 

END 


RCM07800 

RCM07810 

RCM07820 
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APPENDIX B 



PROGRAM RCM 

B. 1 . Program Description 
B.l.l, Computational Grid 

The computational region is divided into 100 cells; the solution grid 
points are odd numbered, e.g., 3, 5, 7 201 with 1 and 203 being tho, 

points where the boundary conditions are specified. The even numbered points, 
2, 4, 6 ..., 202 are intermediate locations where solutions are stored before 
being assigned to the solution grid points. See Fig, (2), 

B.1.2, Data Input 

Data for various ports (exhaust, inlet, etc.) is specified in 
dimensional form in S.I, units (Pascal (N/m^) for pressure, kg/ra^ for 
deru3lty, m/s for velocity etc.). Reference values are also specified in 
like manner. See lines RCM00210 through 00250. 

Initial data is specified through a call to an appropriate 
subroutine, depending on where the calculation is started for a particular 
wave diagram. For the example given in section II on the Spectra Technology 
wave diagram, the computation is started at the point when the high pressure 
gas inlet port just opens. The call for initial data is made to subroutine 
INIT3R, which prescribes data consistent with a solid wall boundary at the 
left and a 'piston* inflow boundary at the right. 

B.1.3, Non ~di mens lonalizat ion 

Non-dimensionalization is carried out in lines 00540 through 00610 
with entropy defined as 

s = Jin 

P ' 
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Note that velocities are all referred to a reference sonic velocity defined 
by 



^ref 



Pref 

^'ref 



B.1.4* Structure 

The main program loop starts at line 00630, for the number of time 
steps specified. The time step is computed according to the appropriate CFL 
condition for the method, and a random number for the time step is generated 
by a call to the function subroutine WDP. 

A secondary loop to define the sequence of local Rieraann problems for 
the time step is set up at line 00750. For each Riemann problem defined, a 
call is made to subroutine GLIMM which i) solves the Riemann problem, and ii) 
samples the solution using the random number generated. The subroutine then 
returns the sampled solution as the parameters PGLIM, RGLIM, UGLIM for the 
pressure, density and velocity respectively. These solutions are initially 
stored in the even numbered intermediate locations on the grid, and are then 
assigned to either the left or the right solution grid point depending on 
whether the random number was In the negative or the positive half of the 
interval respectively. 

A call is then made to one of the modular subroutines structured for 
particular types of wave diagrams, lines 01050-01080, and the others are 
commented out. 

Boundary conditions arci invoked after the call to the modular 
subroutines which return the proper values of the switches SWL and SWR. The 
structure of the boundary condition subroutines is described in section II. 
This sequence completes one pass through the main loop and the process is 
repeated for the number of time steps specified. 
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B . 2 . Example Use of Program RCM 

The program is set up in the following steps: 



i) Line 00150 - 


output device designation. See B.3, 


ii) Line 00190 - 


specify the number of time steps, k , and the 
switches SWL and SWR consistent with where the 
computation is to be started. 


iii) Lines 00210 - 
through 00250 


prescribe flow data for various ports in 
dimensional form. See list of variables for 
explanation of variable names. 


iv) Lines 00490 - 

through 00530 


invoke the proper initial data subroutine and 
comment out the rest. See list of subroutines 
for explanation of subroutine, function subprogram 
names . 


v) Line 00660 


set the interval for number of time steps at which 
a plot of the flow parameters is required. 


vi) Lines 01050 - 

through 01080 


user supplied modular subroutine for particular 
wave diagram to be computed. Comment out the rest 


vii) Line 01190 - 


call to plotting routine should be consistent with 
interval specified in line 0660. 


viii) Lines 02650 - 
through 03700 


identify proper subroutine to prescribe initial 
data (consistent with iv), and specify the data in 
the subroutine in dimensional form. 


ix) Lines 05470 - 

through 05990 


specify plotting parameters, viz., origins of 
plots, scales, number of points to be plotted, 
color of plots, etc. Facility dependent. 




The subroutines PLOTl and PL0T2 given in the 
listing are structured for DISSPLA software 
installed in the facility at NPGS. 


x) Lines 06040 - 

through 07820 


user supplied modular subroutine for wave diagram 
to be computed. 


B.3. Execution 





The program is run in an interactive mode and is invoked through a call 
to DISSPLA, available on most niainframes. After compiling the program 
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(FORTRAN H Extended compiler), the following command executes it: 

DISSPLA filename 

If working at stations equipped with dual screens, the command on line 150 can 
be of the type 

CALL TEK618 Tektronix screen 

If working on a non-graphics terminal, or a single screen station, this 
should be changed to 

CALL COMPRS 

which generates a *DISSPLA METAFILE* to be routed later to either a screen or 
a plotter, e.g., VRSTEC, IBM79, TEK618, etc. Once generated, the metafile can 
be accessed and routed by the command 

DISSPOP device designation 

These are facility dependent commands and should be modified accordingly. 



B . 4 . List of Important Variables (In Alphabetical Order) 



A 

AHEAD 

AMACH 

AL 

AR 

AREF 

ASTAR 

CFLNUM 

DT 

DX 

EPS 

G 

IDIGT 

II 

JCOUNT 

K 

KCOUNT 

N 

N1 



sonic velocity 

sonic speed of head wave of rarefaction fan 
Mach number 

left side sonic speed value for RP 
right side sonic speed value for RP 
reference speed of sound 

speed of sound in ’starred* state of RP solution 
(see Fig. 3) 

CFL number for time step determination 

time step 

grid cell width 

small number for pressure iteration in RP solver 
ratio of specific heats, y 
see WNORM 

argument used in function subprogram PHI equal to either 
0 or 1 
counter 

number of time steps 
counter 

counter for time steps 
counter 
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p 

PA 

PGLIM 

PL 

PR 

PREF 

PSEXIT 

PSINL 

PSINR 

PSOUTn 

PSTAR 

PTOTIN 

QPRINT 

QSTOP 

R 

RA,RB 

RGLIM 

RINL 

RINR 

RL 

RMU 

RR 

RREF 

RTOTIN 

S 

SWL 

SWR 

TCNTCT 

THEAD 

TIME 
TIME REF 
TTOTAL 
U 

UA 

UCNTCT 

UEXMAX 

UGLIM 

UL 

UR 

USTAR 

V^DP 

WL 

WR 

WNORM 

X 

XCNTCT 
XI, XII 



pressure 

flow parameter describing ’a* state in transition functions 

pressure value returned by subroutine GLIMM 

left side pressure value for RP 

right side pressure value for RP 

reference pressure 

static pressure at exit or outlet port 

static pressure for incoming ’piston* flow on left side 
static pressure for incoming ’piston* flow on right side 
n = 1,2,3 - exit static pressures for cycles with more than 
one exhaust port 

pressure in ’starred’ state of RP solution (see Fig. 3) 
total pressure for isentropic inflow 
specification of interval size for output 
maximum number of iterations for solution of Riemann 
problem, (RP) 
density 

flow parameters describing ’a’ and ’b* states in transition 
functions 

density returned by subroutine GLIbIM 

static density for incoming ’piston’ flow on left side 

static density for incoming ’piston’ flow on right side 

left side density for RP 

function of Y 

right side density for RP 

reference density 

total density for isentropic inflow 
entropy 

switch for left boundary 
switch for right boundary 

time taken by contact surface to travel a certain distance 
time taken by head wave of expansion to travel a certain 
distance 

real time in seconds 
reference time 

cumulative non-dimensional time for number of time steps 
velocity 

flow parameter for ’a’ state in transition functions 
velocity of contact surface 

maximum velocity occurring at an outflow boundary 
velocity returned by subroutine GLIMM 
left side velocity for RP 
right side velocity for RP 

velocity in ’starred’ state of RP solution (see Fig. 3) 
value returned by random number generator subprogram 
left shock wave velocity 
right shock wave velocity 

variable used in random number generator subprogram 

space dimension 

location of contact surface 

random numbers scaled to grid cell 
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XREF 

Y 

Z 

ZETA 



- reference length 

- argument used in function subprogram PHI equal to PSTAR 

- argument used in function subprogram PHI equal to either 

PL or PR 

- dummy variable (for initialization purposes in random 

number generator) 



B . 5 . List of Subroutines, Function Subprograms 
B.5.1. Subroutines 



INITl 


prescribes 


initial 


data 


corresponding 


to 


SWL= 1 , 


SWR= 1 ; 




e. g. , shock 


-tube problem 








IN1T2L - 


prescribes 


initial 


data 


corresponding 


to 


SUL=2y 


SWR=1 


INIT2R - 


prescribes 


initial 


data 


corresponding 


to 


SWL= 1 , 


SWR=2 


INIT3L - 


prescribes 


initial 


data 


corresponding 


to 


SWL=3, 


SWR=1 


INIT3R - 


prescribes 


initial 


data 


corresponding 


to 


SWL= 1 , 


SWR=3 


PL0T1,2 - 


graphics subroutines 










GLIMM 


solves the 


Riemann 


problem, samples the 


solution 


[ and 




returns values for 


flow 


parameters 









GE 



DETON 



SPCTRA - 



BCLl 



BCL2 



BCL3 



BCL4 



BCL5 



modular user supplied subroutine to simulate wave diagram 
of General Electric Wave Engine 

modular user supplied subroutine to simulate evacuation of 
detonation chamber 

modular user supplied subroutine to simulate wave diagram 
of Spectra Technology’s Pressure Exchanger 

prescribes boundary conditions (BC’s) corresponding to 
SWL=1, i.e., solid wall on left side 

prescribes BC*s corresponding to SWL=2, i.e«, outflow at 
constant static pressure on left side 

prescribes BC*s corresponding to SWL=3, i.e. , ’piston* 
inflow on left side 

prescribes BC’s corresponding to SVJL=4, i.e., isentropic 
inflow from reservoir on left side 

prescribes BC’s corresponding to SWL=5, i.e., wave ’tuning’ 
on left side 
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BCRl , BCR2, BCR3, BCR4 , BCR5 - prescribe BC’s corresponding to 
SlfR=l ,2,3,4,5 respectively on right side 



B.5.2. Function Subprograms 



PHI(y,z) 

PHIl(PB) 



PSI(PB) 



WDP(II) 



required in iteration procedure for solution of RP 

describes shock transition function, H^a(pb) » 
two states a and b connected by a shock wave (see 
Ref* 6, Ch. Ill) 

describes rarefaction transition function, ^a(pb) > 
for two states a and b connected by a rarefaction wave 
(see Ref. 6, Ch. Ill) 

generates a random number in a van der Corput sequence 
each time it is invoked. Note that it needs to be called 
once from outside the main loop by specifying an argument 
11=1 to initialize IDIGT and l^ORM, returning a value 
of 0 for the dummy variable ZETA, and then a second time 
from within the main loop with an argument 11=0 to return 
a value which is the random number. 
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